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ABSTRACT 

The lack of a consensus bacterial species concept 
greatly hampers our ability to understand and 
organize bacterial diversity. Operational taxonomic 
units (OTUs), which are clustered on the basis of 
DNA sequence identity alone, are the most 
commonly used microbial diversity unit. Although 
it is understood that OTUs can be phylogenetically 
incoherent, the degree and the extent of the phylo- 
genetic inconsistency have not been explicitly 
studied. Here, we tested the phylogenetic signal of 
OTUs in a broad range of bacterial genera from 
various phyla. Strikingly, we found that very few 
OTUs were monophyletic, and many showed 
evidence of multiple independent origins. Using pre- 
viously established bacterial habitats as bench- 
marks, we showed that OTUs frequently spanned 
multiple ecological habitats. We demonstrated that 
ecological heterogeneity within OTUs is caused by 
their phylogenetic inconsistency, and not merely 
due to 'lumping' of taxa resulting from using 
relaxed identity cut-offs. We argue that ecotypes, 
as described by the Stable Ecotype Model, are 
phylogenetically and ecologically more consistent 
than OTUs and therefore could serve as an alterna- 
tive unit for bacterial diversity studies. In addition, 
we introduce QuickES, a new wrapper program for 
the Ecotype Simulation algorithm, which is capable 
of demarcating ecotypes in data sets with tens of 
thousands of sequences. 

INTRODUCTION 

The question of how and whether Ufe organizes itself into 
discrete species units is key to our understanding of how 
diversity originates and is maintained. This question is a 
particularly challenging one for microbiologists because 
unHke plant and animal biologists, we can seldom 



directly observe phenotypic traits that may predict a 
microorganism's ecological niche. Indeed, of the many 
millions of estimated bacterial taxa, the vast majority is 
uncultivable and is known only by DNA sequences (1^). 
This greatly hampers our ability to meaningfully classify 
microbes based on their phenotypes. A key step in 
overcoming this challenge is developing methods to 
organize bacterial DNA sequences into biologically and 
ecologically meaningful taxonomic units. Even whether 
bacterial species exist at all is still a matter of some 
debate among microbiologists (5-8). In the absence of a 
consensus species concept, the most frequently used 
practice for organizing tine scale bacterial diversity is to 
cluster sequences solely on the basis of DNA sequence 
similarity at a conserved locus. Sequence clusters 
dehneated in this manner are termed operational taxo- 
nomic units (OTUs). 

OTUs, clustered based on the 16S rRNA gene, have 
been widely used to approximate bacterial species. 
Although OTUs are expedient for quickly clustering 
large numbers of bacterial sequences, they have several 
significant limitations as a unit of diversity. First, the simi- 
larity cut-off used to define OTUs is arbitrary. Defining 
species by 97% 16S identity is a commonly used rule of 
thumb (9,10), but species so defined are known to encom- 
pass large diversity in genome content, physiology and 
ecology (4,11,12) and are expected to underestimate the 
total diversity present when compared with the accepted 
70% DNA-DNA hybridization threshold (13). Tighter 
thresholds such as 99 or 100% identity have been 
proposed (14) to address this issue, but this only serves 
to highlight the fundamental problem inherent in selecting 
a universal cut-off. It has been shown that regardless of 
the cut-off used, OTUs will not correspond directly to 
existing taxonomic units (15). This is because different 
hneages evolve at different rates; therefore, no universal 
cut-off will capture equivalent units of diversity across aU 
bacterial lineages. 

Another problem with OTUs based strictly on iden- 
tity is that they do not take phylogenetic informa- 
tion into account, as pointed out previously (16-21). 
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Because different lineages evolve at different rates, 
sequence similarity alone is inadequate to infer evolution- 
ary relationships. This problem has been well illustrated in 
the context of gene function prediction based solely on 
sequence similarity (22,23). It is therefore expected that 
similarity-based OTUs wiU contain sequences with mixed 
phylogenetic signal such that an OTU is not guaranteed to 
be monophyletic. However, the extent and degree to which 
this phylogenetic inconsistency exists among OTUs 
remains largely unknown and has not been explicitly 
investigated. Recent studies indicate that bacterial ecolo- 
gical traits in general are phylogenetically conserved 
(17,18,24-31). Assuming phylogenetic niche conservatism, 
we hypothesize that the lack of phylogenetic consistency 
among OTUs wiU be associated with ecological heterogen- 
eity within OTUs. 

These problems with OTUs have led to the development 
of several alternative methods for species demarcation, all 
incorporating phylogeny and having a common ground 
on evolutionary theory. For example, the general mixed 
Yule-coalescent (GMYC) model dehneates the species 
boundary by identifying the transition point from speci- 
ation to coalescent events using a likelihood framework 
(19-21,32). Ecotype Simulation (ES) (33) and AdaptML 
(34) aim to demarcate DNA sequences into ecologically 
cohesive clades (or ecotypes). ES identifies ecotypes by 
comparing the observed pattern of sequence diversity in 
a bacterial community to those of simulated communities 
'evolved' based on the Stable Ecotype Model (8,35). 
AdaptML, by contrast, demarcates ecotypes by inferring 
the evolutionary history of habitat transitions. It identifies 
an ecotype as the largest clade whose members share an 
inferred habitat. 

AdaptML and ES have both been successful in the past 
at predicting bacterial ecotypes from environmental DNA 
sequences. Ecotypes predicted by ES have been confirmed 
as ecologically distinct in isolates from natural 
communities of Bacillus sampled from desert canyons in 
Israel (33) and Death Valley (7), communities of 
Synechococcus from hot springs in Yellowstone 
(11,36,37) and in chnical and environmental isolates of 
Legionella pneumophila (38). AdaptML meanwhile has 
been used to demarcate ecologically distinct clusters in 
marine communities of Vibrio (34) and Desulfobulbus 
(39), as weU as desert soil Bacillus (7). 

By incorporating evolutionary models, GMYC, ES and 
AdaptML carry at least two advantages over OTU clus- 
tering. First, they do not require the selection of an arbi- 
trary sequence identity cut-off. Second, these methods will 
always demarcate species as clades with a single evolution- 
ary origin. 

ES carries the additional advantage of generating pre- 
dictions about the rates of ecotype formation and periodic 
selection within a clade. Periodic selection occurs when an 
individual within an ecotype gains a selective advantage 
within its ecological niche and carries all or nearly all of its 
genotype to fixation within that ecotype (selective sweeps) 
(40,41). This process results in genetic cohesion within 
ecotypes, analogous to the cohesion provided by inter- 
breeding among members of an animal or plant species. 
Ecotype formation occurs when an individual changes its 



ecological niche, thereby releasing itself from the cohesive 
force of periodic selection within its parental ecotype. 
Ecotype formation is therefore analogous to sexual isola- 
tion and speciation among macroorganisms (8,42). 
Estimating the rates of periodic selection and ecotype for- 
mation therefore can shine light on the evolutionary and 
ecological processes that drive microbial diversity. 

Our primary aim for this study was to investigate the 
degree and extent to which OTUs are phylogenetically and 
ecologically inconsistent. We compared and contrasted 
OTUs with ecotypes using both 16S rRNA and protein- 
coding genes. We found surprisingly extensive phylogen- 
etic inconsistency among OTUs, to the extent that only a 
small minority of the OTUs comprised monophyletic 
clades. We also found a large amount of ecological incon- 
sistency among OTUs. Specifically, when tested against 
habitats defined by their ecological parameters, OTUs 
were much more hkely to span multiple habitats and less 
efficient in explaining the ecological variation than were 
ecotypes. 

In addition, we introduce QuickES, a modified version 
of ES that runs much faster than the original version. 
Although an approximation of the complete algorithm, 
this version is capable of generating rough ecotype esti- 
mates for many thousands of sequences, making ecotypes 
a practical alternative unit for microbial diversity studies 
involving large sequence datasets. 



MATERIALS AND METHODS 

Data set, sequence alignment and classification 

The 16SrRNA data set consisted of 116 391 near full- 
length Sanger-sequenced bacterial 16SrRNA sequences 
sampled from 21 different skin sites of 10 human 
subjects (43). The 16SrRNA sequences were aligned 
using the PyNAST algorithm in QIIME (44) and classified 
to the genus level using RDP Classifier, version 2 (45) at 
the default settings. We focused on the 10 most abundant 
genera within the skin data set for subsequent 
analyses. These genera made up ~85% of the 
data set and spanned a broad taxonomic range, including 
the phyla Actinobacteria {Propionibacterium and 
Corynebacterium), Firmicutes {Staphylococcus , 
Streptococcus and Anaerococcus), Bacteroidetes 
(Cloacibacterium) and Proteobacteria (Diaphorobacter, 
Aquabacterium, Acidovorax and Acinetobacter). 

In addition to the 16SrRNA gene, we also analysed 
protein-coding genes including 1025 lisp60 Sanger se- 
quences of the genus Vibrio sampled from a coastal 
marine environment (34) and 132 psaA Sanger sequences 
of the genus Synechococcus sampled from the effluent 
channel of Mushroom Spring in Yellowstone National 
Park (36). The hsp60 encodes a heat shock protein, 
whereas psaA encodes a photosynthetic reaction centre 
protein. Protein-coding sequences were aligned by their 
amino acid sequences using MUSCLE (46) and then con- 
verted back to a DNA alignment using in-house scripts. 
All sequences were retrieved from Genbank, along with 
the sampling data for each sequence. 
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Phylogenetic analysis 

Maximum-likelihood (ML) trees used as input for the ES 
and AdaptML analyses were generated using FastTree 
(47), with the gtr and gamma model engaged. The gtr 
and gamma model was chosen as the best model by 
JModelTest (48). For the assessment of OTU monophyly, 
we generated additional maximum likelihood trees and 
maximum-parsimony trees using RAxML (49), and 
neighbour-joining trees using QuickTree (50). To control 
for uncertainty in the tree topology, we only considered 
clades with >80% bootstrap support in our monophyly 
tests. An OTU was classified as monophyletic if all its 
members shared a single common ancestor to the 
exclusion of other OTUs. Otherwise, it was classified as 
paraphyletic, so long as the last common ancestor of all 
sequences of the OTU had at least 80% bootstrap 
support, and not all the descendants of that ancestral 
node were from the same OTU. We noted, however, 
that paraphyly, so defined, encompassed a broad array 
of phylogenetic patterns. These patterns ranged from 
OTUs that were very close to being monophyletic, 
(i.e. classed as paraphyletic due only to one or two diver- 
gent sequences within an otherwise monophyletic clade), 
to OTUs that were spread across the entire phylogeny. To 
distinguish between these extremes, we computed a 
Taraphyly Index' (PI) (Supplementary Figure SI) for 
each OTU defined by the formula: 

PIOTU = 1 - {NOTU/Nclade) 

Where Nqtu is the number of sequences belonging to 
the OTU, and Ndade is the total number of sequences that 
are descendants of the last common ancestor of the OTU. 
A PI of zero indicates a monophyletic group, whereas a 
PI >0 indicates some degree of paraphyly. 

Although we did not specifically classify any 
OTUs as polyphyletic based on these criteria, we were 
able to determine based on direct observation of 
phytogenies that many OTUs were in fact polyphyletic. 
That is, they were not explainable by a single evolutionary 
origin. 

OTUs 

Unless otherwise specified, OTUs were generated with the 
Uclust de novo clustering algorithm using QIIME (44). 
For the skin 16S rRNA data set, OTUs were generated 
using identity cut-offs of 99.5, 99 and 97%. For the 
protein-coding Vibrio hsp60 and Synechococcus psaA 
data sets, OTUs were generated using cut-offs of 100, 
99, 97, 95, 90 and 85%. It has been shown that different 
OTU clustering algorithms can return very different 
outputs (51). Therefore, to ensure the robustness of our 
analysis of OTU monophyly, additional software 
[MOTHUR (52) and Clusterer (53)] and additional clus- 
tering algorithms (nearest neighbour and average 
neighbour clustering) were also used to cluster the 
sequences. In each case, the default clustering parameters 
were used for aU analyses, except to modify the cut-off 
value, or the clustering algorithm as indicated. 



Demarcation of ecotypes using ES 

We used ES (33) version 0.6 to demarcate ecotypes on two 
abundant genera within the human skin data 
{Aquabacterium and Diaphorobacter), as well as on the 
Vibrio hsp60 and Synechococcus psaA sequences data. 
ES is an algorithm for predicting ecologically homoge- 
neous populations (ecotypes) from sequence data alone, 
without the need for inputting any ecological data, or for 
selecting any similarity cut-off value [see Koeppel et al. 
(33) for a full description of the ES algorithm]. Briefly, 
ES operates by simulating the evolution of a set of 
sequences based on parameterized values for the number 
of ecotypes, as well as the rates of ecotype formation, 
periodic selection and genetic drift. It uses a maximal like- 
hhood framework to estimate values for each of these four 
parameters by fitting the simulated sequences to the 
observed diversity curve (see Supplementary Figure S2 
for examples). Using the best parameter solutions, ES 
then demarcates ecotypes onto a phylogeny generated 
from the same set of sequences, by selecting the most in- 
clusive clades consistent with being a single ecotype. 

The full version of ES is only capable of analysing ~200 
sequences at once within a reasonable time frame. As all of 
these genera contained many more sequences, we used a 
divide-and-conquer approach. Using a guiding tree, we 
subdivided the sequences into clades containing fewer 
than 200 sequences and ran ES separately on each clade. 
We then demarcated ecotypes on the entire tree by finding 
the most inclusive clades consistent with being a single 
ecotype (33). 

QuickES 

To speed up the ES and make it practical to analyse thou- 
sands of sequences, we modified the original ES algorithm 
and created a version called QuickES. QuickES approxi- 
mates the ecotype estimation process but stiU carries the 
key advantages of the standard ES ecotypes over OTUs, in 
that its ecotypes are always monophyletic, and there is no 
need to select an identity cut-off. The improved speed was 
achieved in two primary ways. 

First, we used a divide-and-conquer approach similar to 
the one described earher in the text. In this case, we 
subdivided the sequences into clades such that each 
clade was the most inclusive possible clade in which at 
least 90% of the descendants belonged to the same OTU 
(99% cut-off). Subdividing the data set into clades dra- 
matically improved the speed of analysis because ES scales 
poorly as more sequences are added. 

The second improvement in computation speed came 
from using a rougher estimate of the periodic selection 
and ecotype formation rate parameters. From the set of 
clades obtained in the divide-and-conquer step, we 
selected those clades that contained between 25 and 200 
sequences. Each of these moderately sized clades was then 
analysed using a truncated version of the brute force 
search algorithm from the original ES, to glean a rough 
estimate of the parameter values necessary for ecotype 
demarcation. The very time consuming hill-climbing algo- 
rithm that ES uses to refine the estimates further was 
eliminated. Eliminating this step markedly increases the 
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speed but decreases the reliability of the rate estimates; 
therefore, QuickES should not be used for the final rate 
estimation but only for ecotype demarcations when data 
sets are too large to analyse with the full ES algorithm. 
Having obtained rough estimates of the rate parameters 
for several clades, we then generated global rate estimates 
by computing the mean of the rates from the clades 
analysed. Then, foUowing the basic demarcation 
protocol from the original ES (33), we demarcated 
ecotypes for every clade, by finding the most inclusive 
subclades consistent with being a single ecotype, given 
the global rate estimates. QuickES is freely available 
software and can be downloaded from http:// 
wolbachia.biology.virginia.edu/WuLab/Software.html. 
Detailed instructions of running QuickES are described in 
the Supplementary Methods. 

Demarcation of ecotypes using AdaptML 

We used AdaptML (34) to demarcate ecotypes for the 
marine Vibrio hsp60 data set. Our AdaptML analysis of 
Vibrio returned habitats virtually identical to those of 
Hunt et al. (34), with the exception that we had seven 
habitats instead of six. This is likely due to shght vari- 
ations in tree topology resulting from using different 
tree-building algorithms. We refer to the seventh habitat 
as Hg (Hunt et al. habitats are Ha- Hp). The habitat- 
learning and clustering steps of AdaptML were performed 
using the default settings. 

Benchmarking the performance of OTUs and ecotypes in 
explaining the ecological variance in the vibrio data set 

We used the methods described in (19) to compare OTUs 
and ecotypes in their ability to account for the ecological 
variation in the Vibrio data set. Vibrio sequences were 
associated with two ecological measurements: the size of 
the particle from which the Vibrio sequences were isolated 
(<1 |Am, 1-5 |im, 5-63 |im, >63 |im) and the season when 
the samples were collected (spring and fall) (34). We first 
transformed the particle size and season categories into 
quantitative variables using multiple correspondence 
analysis implemented in the 'ade4' package (54). Then 
we carried out redundancy analysis implemented in the 
'vegan' package in R to estimate the amount of variation 
in the ecological parameters that could be explained when 
sequences were grouped by either OTUs or ecotypes. 
Statistical significance was assessed using permutation 
tests. Akaike information criterion was used to evaluate 
the performance of the different species dehneation 
models. 

Estimation of periodic selection and ecotype 
formation rates 

The complete ES algorithm generates estimates of the 
rates of periodic selection and ecotype formation. We per- 
formed ES analyses individually on all major subclades of 
the genera Aquabacterium, Diaphorobacter and Vibrio and 
tested whether the mean rates of ecotype formation and 
periodic selection were different between genera using 
Student's /-tests. Rates within the same genus were 
grouped into high, medium and low categories using the 



Tukey-Kramer test as described in Supplementary 
Methods. ES estimated rates in units of events per nucleo- 
tide substitution. Rates were log-transformed before stat- 
istical analysis so as to more closely approximate a normal 
distribution. 

RESULTS 

Extensive microdiversity within clades 

Consistent with our expectations, and with findings 
observed in other microbial habitats, we observed exten- 
sive microdiversity in the data sets we analysed 
(Supplementary Figure S2). The number of 16SrRNA 
OTUs in each genus showed a dramatic flare-up between 
the 98 and 99% cut-off levels (Supplementary Figure 
S2A-E). The presence of such a 'hockey-stick' pattern 
has previously been observed in natural bacterial popula- 
tions (55,56) and is considered typical of human 
associated microbial populations (57,58). This pattern is 
consistent with the Stable Ecotype Model (8,33), which 
predicts that ephemeral microdiversity should be present 
within bacterial communities, as bacteria undergo neutral 
divergence between periodic selective sweeps. 

The protein-coding sequences of the Vibrio and 
Synechococcus data sets displayed similar flare-ups, 
though at shghtly lower sequence identity thresholds 
(~97% in both cases) (Supplementary Figure S2 Fand 
G). This reflects the more rapid evolution of the protein- 
coding hsp60 and psaA genes compared with the 
16SrRNA gene. 

Extensive and pronounced paraphyly and polyphyly 
among OTUs 

To deal with the potential uncertainty in the tree topology, 
we only considered clades that were weU supported (boot- 
strap values >80) in our monophyly analyses. Our analysis 
of the predominant skin bacterial genera revealed that 
strikingly few OTUs are monophyletic (Table 1). At the 
97% identity level, no >75% of the OTUs were monophy- 
letic groups in any of the genera analysed. At the 99% 
identity level, fewer than 60% of OTUs in each genus were 
monophyletic. The percentages were far smaller in most 
genera. The results were even more remarkable among 
larger OTUs (those containing at least 50 sequences): 
<67% of the large OTUs in each genus were 
monophyletic at the 97% cut-off. At the 99% cut-off, 
fewer than 25% of the large OTUs in each genus were 
monophyletic groups. In fact, in five of the 10 genera 
analysed {Acidovorax, Acinetobacter, Aquabacterium, 
Cloacibacterium and Diaphorobacter), none of the large 
99% OTUs were monophyletic. 

To measure the degree of paraphyly among these 
OTUs, we computed the PI (see 'Materials and 
Methods' section for details) for all of the 99% OTUs in 
each genus. A PI of 0 indicates a monophyletic group, 
whereas a PI close to 1 indicates substantial paraphyly. 
Surprisingly, large numbers of OTUs of all sizes across 
all genera had mixed phylogenetic signal, some extensively 
(i.e. their PI was close to 1.0, Figure 1). We observed 
similar patterns among OTUs clustered at 97 and 99.5% 
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Table 1. Phylogenetic heterogeneity among 16S rRNA OTUs of skin data set 



Genus OTU identity OTUs >1 Sequence OTUs >50 Sequences 

thresliold 

Number % Monophyletic Number % Monophyletic 

of OTUs of OTUs 



Acidovorax 


97% 


2 


0.00% 


1 


0.00% 




99% 


3 


0.00% 


1 


0.00% 




99.5% 


4 


0.00% 


2 


0.00% 


Acinetohacter 


97% 


8 


12.50% 


3 


0.00% 




99% 


31 


25.81% 


8 


0.00% 




99.5% 


45 


31.11% 


1 


0.00% 


Anaerococcus 


97% 


31 


51.61% 


4 


0.00% 




99% 


64 


37.50% 


4 


25.00% 




99.5% 


71 


11.27% 


3 


0.00% 


Aquahacterium 


97% 


5 


60.00% 


1 


0.00% 




99% 


5 


60.00% 


1 


0.00% 




99.5% 


12 


41.67% 


1 


0.00% 


Cloacihacterium 


97% 


2 


50.00% 


1 


0.00% 




99% 


7 


57.14% 


2 


0.00% 




99.5% 


44 


11.36% 


3 


0.00% 


Corvnehacterium 


97% 


53 


33.96% 


22 


13.64% 




99% 


160 


25.00% 


44 


2.27% 




99.5% 


323 


18.89% 


45 


0.00% 


Diaphorohacter 


97% 


1 


0.00% 


1 


0.00% 




99% 


1 


0.00% 


1 


0.00% 




99.5% 


11 


18.18% 


1 


0.00% 


Propionihacterium 


97% 


3 


66.67% 


3 


66.67% 




99% 


7 


42.86% 


4 


25.00% 




99.5% 


88 


23.86% 


3 


0.00% 


Staphylococcus 


97% 


6 


16.67% 


3 


0.00% 




99% 


59 


18.64% 


16 


12.50% 




99.5% 


230 


24.78% 


12 


0.00% 


Streptococcus 


97% 


12 


75.00% 


5 


60.00% 




99% 


50 


22.00% 


7 


14.29% 




99.5% 


127 


14.96% 


4 


0.00% 



This table displays the number of monophyletic OTUs in each genus at three different identity thresholds (97, 99 and 99.5%). Only 
OTUs containing at least two sequences and meeting the support criteria were considered, as a single sequence is monophyletic by 
definition. The effect was more pronounced among larger OTUs (OTUs containing at least 50 sequences, right-hand columns). 



thresholds (Suppleinentary Figure S3). This resuh suggests 
that phylogenetic incoherence among OTUs is far more 
pronounced and pervasive than is generally recognized. 

Given the extreme deviation froin monophyly that we 
observed with the PI scores, we mapped several OTUs 
onto phylogenies to observe their phylogenetic patterns 
directly. We observed that many of the OTUs appear to 
be polyphyletic because they required multiple independ- 
ent evolutionary origins to explain their distribution 
across the phylogeny (Figure 2). Extensive paraphyly 
was observed, regardless of the phylogenetic or OTU clus- 
tering methods used (Table 2). Putative ecotypes, by 
contrast, always inapped to monophyletic clades 
(Figure 2) because phylogenetic information was taken 
into account during the demarcation process. 

We next checked whether the phylogenetic heterogen- 
eity we had observed among OTUs based on 16S rRNA 
also existed in OTUs of protein-coding genes. OTUs in the 
marine Vibrio data set, clustered based on similarity at the 
hsp60 locus, showed a similar, if less extreine phylogenetic 
pattern. There was considerably more monophyly among 
OTUs in this data set than in the skin data set. Over 80% 
of the OTUs at all cut-offs except 99% were monophyletic 
(Supplementary Table SI). At the 99% cut-off, however, 
only 74.56% of the OTUs were monophyletic groups. 



Although the inajority of the non-monophyletic OTUs 
were paraphyletic (i.e. explainable by a single evolutionary 
origin), several OTUs did appear to be polyphyletic 
(Figure 3A) (see 97% OTU-6 and OTU-7). The 
Yellowstone Synecliococcus psaA sequences clustered 
into only three OTUs (97% cut-off). Two of the three 
were monophyletic, but the third and largest was paraphy- 
letic (Figure 3B). 

Extensive ecological heterogeneity among OTUs 

It has been estabhshed that OTUs contain a great deal of 
ecological diversity (11). However, this is usually assumed 
to be a result of defining OTUs too broadly. Under this 
hne of reasoning, it inight be assumed that simply narrow- 
ing the breadth of OTUs, either by using a more stringent 
identity cut-off or by using a marker less conserved 
than 16S rRNA, can address this problem in microbial 
ecology studies. We have undercut this reasoning by 
deinonstrating that OTUs are phylogenetically inconsist- 
ent. Our results suggest that the problein runs deeper than 
simply using a cut-off that is too relaxed. We hypothesized 
that the lack of phylogenetic consistency among OTUs 
will be associated with ecological heterogeneity within 
OTUs, even when a more rapidly evolving protein- 
coding gene was used to cluster them. 
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Figure 1. OTU paraphyly is pervasive and pronounced. This graph plots OTU size against PI for all 99% 16S rRNA OTUs among 10 genera. PI 
values of 0.0 indicate monophyletic groups, whereas a PI close to 1 indicates substantial paraphyly. Genus classifications of OTUs are colour coded 
as indicated in the key. 



We used the previously established marine Vibrio (34) 
and hot spring Synechococcus (36) AdaptML habitat types 
as benchmarks for determining ecological homogeneity. 
We reasoned that populations belonging to a single 
habitat could be considered more ecologically homoge- 
neous than populations spanning multiple habitats. We 
observed that in the Vibrio data, OTUs (97% identity 
cut-off at the hsp60 locus) frequently spanned several of 
the habitats predicted by AdaptML (Figure 3A). 

For the entire Vibrio data set, we found that only 55.1 % 
of the OTUs belonged to just one AdaptML habitat. In 
comparison, 77.2% of ecotypes generated by ES were 
either identical to or nested within the AdaptML 
ecotypes, as expected (35) (and see 'Discussion' section). 
A Student's ?-test confirmed that the mean number of 
habitats spanned by Vibrio OTUs was significantly 
greater than the mean number spanned by Vibrio ES 
ecotypes {P < 0.0002). This indicated that ES ecotypes 
showed greater ecological homogeneity than the 97% 
OTUs in the Vibrio data set. 

The Synechococcus analysis showed a similar pattern. 
Our ES analysis predicted 12 non-singleton ES ecotypes 
in the Synechococcus data set (Figure 3B), which appear to 
be ecologically distinct based on the temperature of the 
sample site (36). Representatives of several ecotypes were 
also found to be predominant at different depths of the 
microbial mat, suggesting potential ecological distinctness 
based on light and O2 concentrations (36). Each of the ES 
ecotypes corresponded to a single AdaptML habitat, 
though, as expected, some AdaptML ecotypes contained 
multiple ES ecotypes (Figure 3B). In comparison, the 
largest OTU spanned all three AdaptML habitats. 

The size of an OTU is expected to have a large impact 
on its level of ecologically heterogeneity. For example, for 



the same data set, 97% OTUs will be larger than 99% 
OTUs (i.e. containing more sequences) and thus ecologic- 
aUy will be more heterogeneous. To control for the poten- 
tial size difference between OTUs and ecotypes, we 
foUowed the approach of Powell et al. (19) and bench- 
marked the performance of OTUs of different cut-offs 
and ecotypes in explaining the ecological variance in the 
Vibrio data set. Our results show that ecotypes (ES or 
AdaptML) significantly outperformed OTUs in account- 
ing for variation in the ecological parameters (samphng 
particle size and season) (Table 3). This was the case, re- 
gardless of the size of the OTUs (at 97, 99 or 100% 
identity cut-offs). Although 100% OTUs explained the 
most variance, as expected, it came with a high cost 
associated with the large number of classes. Ecotypes 
most efficiently explained the variance according to the 
Akaike information criterion (8A1C =181 for AdaptML 
and 5AIC = 49 for ES against 100% OTUs, respectively). 
It is not surprising that AdaptML performed better than 
ES in this test because AdaptML uses both the sequences 
and the ecological information to demarcate ecotypes, 
whereas ES only uses sequences. Our results are in 
agreement with the Powell et al. study (19) showing that 
evolutionary theory-based approaches outperform oper- 
ational approaches in producing ecologically meaningful 
diversity units. 

Universal identity cut-offs fail to capture all 
putative ecotypes 

Having shown that 97% cut-off OTUs were ecologically 
heterogeneous, we next sought to determine whether any 
single identity cut-off could have generated OTUs similar 
to the putative ecotypes designated by the ES algorithm. 
We binned the ecotypes by their minimum pairwise 
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Figure 2. Extensive paraphyly and polyphyly among OTUs. Maximum likelihood trees of representative subclades of the genera (A) Aquabacterium 
and (B) Diaphorohacter. OTU generated using the 99% identity cut-off are shown with the putative ecotypes (PE) demarcated by ES. Internal nodes 
with >80% bootstrap support are highlighted with red circles. 



5182 Nucleic Acids Research, 2013, Vol. 41, No. 10 



Table 2. Phylogenetic heterogeneity of OTUs is robust to methodology 
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This table displays the percentage of monophyletic 16S rRNA OTUs in the genus Aquahacterium at three different identity thresholds (97, 99 and 
99.5%). As in Table 1, only OTUs containing at least two sequences and meeting support criteria were considered in the percentage computations, 
though the total number of OTUs is displayed in the n column. Different phylogenetic methods (columns) and different OTU clustering algorithms 
(rows) were tested. Cells display the percentage of OTUs that were monophyletic clades. 
ML, Maximum likelihood; MP, Maximum parsimony; NJ, Neighbour joining. 



sequence identity and plotted the number of ecotypes in 
each sequence identity bin for the Aquahacterium and 
Vibrio genera (Figure 4). Figure 4 shows that ecotypes 
display a wide range of sequence identities, and no 
universal OTU cut-off can be applied to capture all 
these ecotypes. This held true, regardless of whether the 
complete ES method or QuickES was used. For example, 
in the Aquahacterium genus, a large majority (94.2%) of 
the ecotypes had a minimum pairwise sequence identity 
>99% (Figure 4 A left panel). This means that if these 
sequences were clustered into OTUs using a 99% 
identity cut-off, many of those OTUs might contain 
sequences from multiple ecotypes. At the other end of 
the spectrum are the 14 ecotypes (5.8%) whose 
minimum sequence identity is <99%. A 99% identity 
OTU cut-off would subdivide these ecotypes into 
multiple OTUs. 

Rates of periodic selection vary within and between genera 

One advantage of ES over other methods is that it can be 
used to gain insight on the processes of microbial 
diversification. We used the complete ES to compare the 
rates of periodic selection and ecotype formation within 
and between genera. We generated estimates of the 
rates of periodic selection and ecotype formation for 
subclades within each genus. Consistent with the Stable 
Ecotype Model, we found that, in each case, the rate 
of periodic selection was estimated to be greater than 
the rate of ecotype formation. In Aquahacterium and 
Diaphorohacter, the median rate of periodic selection 



among subclades was around twice the median rate of 
ecotype formation (2.04 times and 1.50 times greater, 
respectively). Interestingly, in Vihrio, the median rate of 
periodic selection was 26.5 times higher than the median 
rate of ecotype formation (Table 4). Our analysis 
revealed no statistically significant differences between 
the Vihrio, Aquahacterium and Diaphorohacter genera in 
their rates of ecotype formation (Student's /-test: 
Vihrio-Diaphorohacter. P < 0.298; Vihrio- Aquahacterium: 
P < 0.097; Aquabacterium-Diaphorohacter: P < 0.997). 
However, we did detect significant differences in 
the rate of periodic selection between genera. Specifically, 
both Vihrio and Aquahacterium showed a significantly 
higher rate of periodic selection than Diaphoro- 
hacter (Student's ?-test: /'< 0.0059 and /"< 0.0139, 
respectively). 

There were surprisingly large variations in the rates of 
periodic selection and ecotype formation (Supplementary 
Figures S4 and S5) within each genus. For example, both 
Aquahacterium and Vihrio contained clades with greatly 
elevated periodic selection rates. In Aquahacterium, the 
rate estimated for one clade was nearly 8-fold higher 
than that of the other clades (Supplementary Figure 
S4B). In Vihrio, two clades showed rates 16-fold and 
160-fold higher above the average, respectively 
(Supplementary Figure S4C). The rate of ecotype forma- 
tion within Vihrio was also highly variable, with one group 
of clades showing ecotype formation rates as much as 
10-fold higher than the rest of the genus (Supplementary 
Figure S5C). 
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Figure 3. OTUs and Ecotypes show distinct habitat associations. An ML tree of a subset of the Vibrio h.sp60 sequences (A) and a neighbour-joining 
tree of the full set of Synechococcus psaA sequences (B). OTUs, ES ecotypes and AdaptML ecotypes are shown. Note that the formatting in the OTU 
column and the AdaptML column is different. In the OTU column, all leaves marked by the same color belong to the same OTU. In the AdaptML 
column, different colours denoted different habitats. Each distinct colour bar is its own ecotype, whereas bars of the same colour are ecotypes co- 
occurring in the same habitat. 
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Validation of QuickES 

We tested QuickES on a three-gene concatenation of se- 
quences from soil Bacillus simplex and Bacillus subtilisj 
licheniformis. These sequences were isolated from the 
Negev Desert in Israel and had previously been analysed 
using the complete ES algorithm (33). For each data set, 
we ran three QuickES trials, selecting the results from the 
trial whose parameter solutions gave the best hkelihood 
value (best of three). We then repeated this procedure 



Table 3. Performance of OTUs and ecotypes in explaining ecological 
variation in the Vibrio data set 



Model 


Number 


Variance 


AIC 




of classes 


explained 




OTU 97% 


68 


42%* 


764 


OTU 99% 


187 


60%* 


649 


OTU 100% 


382 


73%* 


639 


ES 


190 


62%* 


590 


QuickES 


156 


58%* 


609 


AdaptML 


214 


69%* 


458 



The variance explained was calculated by dividing the constrained 
inertia by the total inertia. The two methods that returned the lowest 
AIC scores are highlighted in bold. Asterisk denotes the significance of 
P < 0.005 by permutation tests. 
AIC, Akaike information criterion. 



twice, resulting in a total of three repHcate runs. 
Supplementary Figures S6 and S7 show that ecotypes 
demarcated in three QuickES runs were very similar to 
each other and also to those predicted by the fuU analysis. 

We also benchmarked the performance of QuickES 
ecotypes in its ability to explain the ecological variance 
in the Vibrio dataset (Table 3). Grouping sequences by 
QuickES ecotypes accounted for 58% of ecological vari- 
ation (P< 0.005). QuickES ecotypes were much better at 
explaining the ecological variance than all the OTUs we 
tested (97%, 99% and 100%, 8AIC>30), although they 
performed slightly worse than the ecotypes demarcated 
with the complete ES (6AIC = 19). 

DISCUSSION 

Although great strides have been made in our understand- 
ing of bacterial diversity in recent years, the field has been 
challenged by the lack of an agreed on unit of diversity 
analogous to the biological species of animals and plants. 
Although identity-based OTUs provide a 'quick and dirty' 
approach to quantifying bacterial diversity and are very 
useful, they are no substitute for coherent and meaningful 
units of bacterial ecology and evolution (8,59,60). 

Our results illustrate several problems with using OTUs 
to approximate bacterial species. Although these problems 
are generally known, we have shown here that they are far 
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Figure 4. Putative ecotypes are not captured by any single sequence identity cut-off Graphs display the number of ecotypes whose minimum 
pairwise sequence identity fall into each of the displayed bins in the skin Aquabacterium (A) and the marine Vibrio (B) data sets. 
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Table 4. Ecotype formation and periodic selection rates of three genera 



Genus n Ecotype formation rate Periodic selection rate 



Mean Standard Median Mean Standard Median 

deviation deviation 



Aquabacterium 20 0.179 0.092 0.458 6.969 18.950 0.934 

Diaphorohacter 19 0.166 0.130 0.340 0.707 0.472 0.510 

Vibrio 10 0.312 0.370 0.105 88.467 239.166 2.780 



ES (original version) estimates of the mean and median periodic selection and ecotype formation rates for skin Aquabacterium and 
Diaphorohacter, and for marine Vibrio. The n values are the number of subclades for which the rates were independently calculated. 



more serious than has previously been appreciated. First, 
paraphyly and polyphyly among OTUs are far more 
extensive and pronounced than expected. We observed 
this pattern in a wide variety of genera spanning 
muhiple bacterial phyla, indicating that this problem is 
widespread and not taxon specific. Our results showing 
significant numbers of paraphyletic and polyphyletic 
OTUs at various sequence identity thresholds further re- 
inforce the point that simply narrowing the identity cut- 
off will not correct the phylogenetic inconsistency among 
OTUs. Although the acceptability of paraphyletic taxa in 
systematics is debatable (8,61,62), polyphyletic groups are 
generally considered to be unacceptable as true taxa. This 
is because polyphyletic groups by definition lack a single 
evolutionary origin. Most modern species concepts require 
individuals within the same species to share a single 
evolutionary lineage (63,64) and therefore a single evolu- 
tionary origin. 

Second, we show that extensive phylogenetic inconsist- 
ency is associated with extensive ecological heterogeneity 
in OTUs. This is distinct from ecological heterogeneity 
caused by using OTUs that are too broadly defined. Just 
as phylogenetic inconsistency will persist, regardless of 
identity cut-off, simply narrowing the identity cut-off 
might not produce ecologically coherent OTUs. We 
demonstrated that ecotypes outperform OTUs in explain- 
ing the ecological variance in the sequence data. This 
finding is particularly relevant to human microbiome 
research, in which 16S rRNA OTUs are commonly used 
to compare the composition of microbial communities 
between healthy and diseased individuals (65). The 16S 
rRNA OTUs at 97 (66,67), 98 (58) and 99% (43,68) 
identity cut-offs have been used as units of microbial di- 
versity. Our results suggest that such OTUs may not be 
ecologically homogeneous and therefore can be problem- 
atic for association studies. Specifically, the use of eco- 
logically heterogeneous units could add noise when 
investigating the associations between OTUs and health 
states, resulting in false negatives. 

An additional problem with OTUs is that the identity 
cut-off is arbitrarily defined and subjectively appUed. The 
threshold selected can greatly affect estimates of both the 
numbers and composition of species in a community. 
Ecotypes, as predicted by either ES or AdaptML, 
require no arbitrary cut-off, and therefore can be 
compared consistently across taxa. We demonstrated 
that ecotypes can be variable in the amount of sequence 



diversity they contain such that no single universal identity 
threshold will accurately capture the ecotypes in a com- 
munity. We calculated alpha diversity (diversity within 
habitats) and beta diversity (diversity among habitats) 
indices for the skin data set using both OTUs and 
ecotypes (data not shown). The overall alpha diversity 
values measured using ecotypes are higher than those 
measured using 16SrRNA OTUs (97 or 99% cut-off), 
suggesting that the routinely used 16SrRNA units might 
underestimate the diversity. However, we noticed no 
apparent or significant differences between beta diversity 
measures derived from ecotypes and OTUs. 

One benefit of the complete ES algorithm over other 
methods is that it provides estimates of the rates of 
ecotype formation and periodic selection. These rate esti- 
mates can be used to gain insight into the mode and tempo 
of bacterial diversification. Periodic selection events purge 
the genetic diversity within a population and lead to adap- 
tation within a species lineage (anagenesis) (40,41). 
Ecotype formation (cladogenesis) occurs when an individ- 
ual becomes sufficiently ecologically distinct that it is no 
longer vulnerable to periodic selection events occurring in 
its parent population. Determining the relative rates of 
periodic selection and ecotype formation can help us dis- 
tinguish between different models of bacterial speciation 
and evolution (8,69). For example, the Stable Ecotype 
Model proposes recurring periodic selective sweeps 
between rare ecotype formation events, and therefore 
predicts that anagenesis is the most dominant mode of 
adaptive evolution (8). In contrast, the Species-Less 
Model features a high rate of species turnover, with 
frequent cladogenesis and almost no anagenesis, because 
each species is hkely to go extinct before its first periodic 
selection event (69). Our ES analyses indicate that periodic 
selection happens much more frequently than ecotype for- 
mation, suggesting that at least in the bacterial lineages 
analysed in our study, anagenesis is the dominant mode of 
adaptive evolution. 

For a lineage that is recently formed, either as a result of 
invading a new habitat or via horizontal transfer of a 
niche-expanding gene (35), we might expect an elevated 
rate of anagenesis as the lineage adapts to the conditions 
of its new environment. Highly elevated periodic selection 
rates Uke those we observed in some subclades of Vibrio 
and Aquabacterium therefore might be an indication of 
their recent expansion into a novel ecological niche. This 
is consistent with a punctuated equilibrium mode of 
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evolution, in which evolutionary change occurs in rare but 
rapid bursts following major changes in a hneage's envir- 
onment (70). This hypothesis could be tested using com- 
parative genomic and ecological analyses between sister 
clades with high and low periodic selection rates, for 
example, to look for genes that are undergoing positive 
selection. 

Capable of processing only hundreds of sequence at a 
time, the original version of ES would be unable to prac- 
tically analyse massive next generation sequencing data 
sets. We have demonstrated here that QuickES, by using 
a divide and conquer approach coupled with parameter 
approximation, can demarcate ecotypes in datasets with 
tens of thousands of sequences. QuickES generates only 
very rough estimates of the periodic selection and ecotype 
formation rates, and we do not recommend QuickES be 
used to draw conclusions of this type. In addition to the 
QuickES package introduced here, a new version of 
the ES algorithm is under development that should 
allow the fuU ES algorithm to directly analyse many thou- 
sands of sequences at once (Frederick M. Cohan, Danny 
Krizanc, personal communication). AdaptML can already 
analyse thousands of sequences simultaneously but 
requires ecological data to estimate ecotypes. ES needs 
only sequence data to predict ecotypes and therefore is 
advantageous when little or no ecological data is available. 
When possible, the best practice is to use both AdaptML 
and ES together to take advantages of their complemen- 
tary benefits, with AdaptML demarcating putative 
ecotypes based on known ecological parameters, and ES 
then subdividing them based on evolutionary models, as 
we and others have demonstrated (7,36,37). Because 
ecotypes incorporate evolutionary and ecological models, 
they are evolutionarily and ecologically more consistent 
than OTUs. Given the numerous advantages of ecotypes 
over OTUs, we advocate for using ecotype as an alterna- 
tive unit of microbial diversity. 
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